#Multivariate stats of output from skyline daily - all proteins and all oysters
setwd('/Users/emmatimminsschiffman/Documents/Dissertation/proteomics/DB post-genome/OT expression comparison/Skyline daily/Re-analysis 072613')

sky.dat<-read.csv('skyline for NMDS.csv', header=T, row.names=1)
sky.dat2<-read.csv('skyline for NMDS pCO2 only.csv', header=T, row.names=1)

setwd('/Users/emmatimminsschiffman/Documents/Dissertation/proteomics/DB post-genome/OT expression comparison')
treatment<-read.csv('treatment groupings.csv', header=T, row.names=1)

sky.t<-t(sky.dat)
sky.t2<-t(sky.dat2)

#there are 2 0's in the dataset so need to log(x+1) transform
skydat.tra<-(sky.t+1)
skydat.tra<-data.trans(skydat.tra, method='log', plot=F)

#NMDS
skydat.nmds<-metaMDS(skydat.tra, distance='bray', k=2, trymax=100, autotransform=F)
skydat.vec<-envfit(skydat.nmds$points, skydat.tra, perm=100)
ordiplot(skydat.nmds, choices=c(1,2), type='text', display='sites', xlab='Axis 1', ylab='Axis 2')
plot(skydat.vec, p.max=0.01, col='blue', cex=0.5)
ordihull(skydat.nmds, treatment[,3], label=F, col='blue')

skydat.high<-metaMDS(skydat.tra[1:8,], distance='bray', k=2, trymax=100, autotransform=F)
high.vec<-envfit(skydat.high$points, skydat.tra[1:8,], perm=100)
fig.high<-ordiplot(skydat.high, type='none', ylim=c(-0.02, 0.02), xlim=c(-0.02, 0.03), main='C')
points(fig.high, 'sites', col='orangered2', pch=c(rep(17,4), rep(2,4)))
ordihull(skydat.high, treatment[1:8,2], label=F, col='grey47')
legend(0.009, 0.02, legend=c('unstressed', 'stressed'), col='orangered2', pch=c(17,2))

skydat.low<-metaMDS(skydat.tra[9:16,], distance='bray', k=2, trymax=100, autotransform=F)
low.vec<-envfit(skydat.low$points, skydat.tra[9:16,], perm=100)
fig.low<-ordiplot(skydat.low, type='none', main='B', xlim=c(-0.01, 0.02), ylim=c(-0.01, 0.015))
points(fig.low, 'sites', col='dodgerblue2', pch=c(rep(15,4), rep(0,4)))
ordihull(skydat.low, treatment[1:8,2], label=F, col='grey47')
legend(0.0075, 0.015, legend=c('unstressed', 'stressed'), col='dodgerblue2', pch=c(15,0))

skydat2.tra<-(sky.t2+1)
skydat2.tra<-data.trans(skydat2.tra, method='log', plot=F)
skydat.co2<-metaMDS(skydat2.tra, distance='bray', k=2, trymax=100, autotransform=F)
co2.vec<-envfit(skydat.co2$points, skydat2.tra, perm=100)
fig.co2<-ordiplot(skydat.co2, type='none', main='A', xlim=c(-0.01, 0.015), ylim=c(-0.01, 0.01))
points(fig.co2, 'sites', col=c(rep('orangered2', 4), rep('dodgerblue2', 4)), pch=c(rep(17,4), rep(15,4)))
ordihull(skydat.co2, treatment[5:12,1], label=F, col='grey47')
legend(0.00505, 0.01, legend=c('400 µatm', '2800 µatm'), col=c('dodgerblue2', 'orangered2'), pch=c(15,17))

plot.nmds<-par(mfrow=c(2,2))

#anosim - low vs. high pCO2
#normalize dataset by row
co2.row<-data.stand(sky.t2, method='total', margin='row', plot=F)
#create dissimilarity matrix using Bray-Curtis coefficient
co2.d<-vegdist(co2.row, 'bray')

#anosim for pCO2 grouping
co2.anos<-anosim(co2.d, grouping=treatment[5:12,1])
summary(co2.anos)
#R=-0.1875, p=0.123

#anosim - MS low pCO2
dat.row1<-data.stand(sky.t[9:16,], method='total', margin='row', plot=F)
dat.d1<-vegdist(dat.row1, 'bray')
low.anos<-anosim(dat.d1, grouping=treatment[9:16,3])
summary(low.anos)
#R=0.1042, p=0.326
#Error in out[i, ] <- tmp[[i]] : 
  #number of items to replace is not a multiple of replacement length

#anosim - MS high pCO2
dat.row2<-data.stand(sky.t[1:8,], method='total', margin='row', plot=F)
dat.d2<-vegdist(dat.row2, 'bray')
high.anos<-anosim(dat.d2, grouping=treatment[1:8,3])
summary(high.anos)
#R=0.08333, p=0.302
#same error

#anosim - all
dat.row<-data.stand(sky.t, method='total', margin='row', plot=F)
dat.d<-vegdist(dat.row, 'bray')
all.anos<-anosim(dat.d, grouping=treatment[,3])
summary(all.anos)
#R=0.02865, p=0.352


#ANOSIM for highly signficant proteins based on eigenvector loadings
lowms.sig<-read.csv('sig proteins 400 mechstress.csv', header=T, row.names=1)
lowms.t<-t(lowms.sig)
lowms.row<-data.stand(lowms.t, method='total', margin='row', plot=F)
lowms.d<-vegdist(lowms.row, 'bray')
anos.lowms<-anosim(lowms.d, grouping=treatment1)
summary(anos.lowms)
#R=0.09375, p=0.235

highms<-read.csv('sig proteins 2800 mechstress.csv', header=T, row.names=1)
highms.t<-t(highms)
highms.row<-data.stand(highms.t, method='total', margin='row', plot=F)
highms.d<-vegdist(highms.row, 'bray')
anos.highms<-anosim(highms.d, grouping=treatment1)
summary(anos.highms)
#R=0.1354, p=0.247

pco2<-read.csv('sig proteins OA.csv', header=T, row.names=1)
pco2.t<-t(pco2)
pco2.row<-data.stand(pco2.t, method='total', margin='row', plot=F)
pco2.d<-vegdist(pco2.row, 'bray')
anos.pco2<-anosim(pco2.d, grouping=treatment2)
summary(anos.pco2)
#R=0.5208, p=0.022

treatment1<-c(rep('NS',4), rep('MS',4))
treatment2<-c(rep('High', 4), rep("Low", 4))

#heat maps for highly significant proteins based on eigenvector loadings

OAprot.names<-read.csv('sig OA NMDS proteins with gene names.csv', header=T, row.names=1)
OAprot.log<-data.trans(OAprot.names, method='log', plot=F)
pheatmap(OAprot.log, cluster_rows=T, cluster_cols=T, clustering_distance_rows='euclidean', clustering_distance_cols='euclidean', clustering_method='average', show_rownames=F, show_colnames=F)
#2800 oysters are 4 on the left, 400 oysters are 4 on right